//lc: size of mesh
lc = 9000;

//points
Point(1) = {-50000, -50000, 30000, lc};
Point(2) = {50000, -50000, 30000, lc};
Point(3) = {50000, 50000, 30000, lc};
Point(4) = {-50000, 50000, 30000, lc};
Point(5) = {-50000, -50000, 0, lc};
Point(6) = {50000, -50000, 0, lc};
Point(7) = {50000, 50000, 0, lc};
Point(8) = {-50000, 50000, 0, lc};
Point(9) = {-50000, -50000, -50000, lc};
Point(10) = {50000, -50000, -50000, lc};
Point(11) = {50000, 50000, -50000, lc};
Point(12) = {-50000, 50000, -50000, lc};

Point(13) = {-500, -1000, -250, lc/200};
Point(14) = {500, -1000, -250, lc/200};
Point(15) = {500, 1000, -250, lc/200};
Point(16) = {-500, 1000, -250, lc/200};
Point(17) = {-500, -1000, -2250, lc/200};
Point(18) = {500, -1000, -2250, lc/200};
Point(19) = {500, 1000, -2250, lc/200};
Point(20) = {-500, 1000, -2250, lc/200};

//测线点

Point(21) = {0,-1000,0,lc/300};
Point(22) = {0,1000,0,lc/300};




//lines
Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {3,4};
Line(4) = {4,1};
Line(5) = {5,6};
Line(6) = {6,7};
Line(7) = {7,8};
Line(8) = {8,5};
Line(9) = {9,10};
Line(10) = {10,11};
Line(11) = {11,12};
Line(12) = {12,9};
Line(13) = {1,5};
Line(14) = {2,6};
Line(15) = {3,7};
Line(16) = {4,8};
Line(17) = {5,9};
Line(18) = {6,10};
Line(19) = {7,11};
Line(20) = {8,12};

Line(21) = {13,14};
Line(22) = {14,15};
Line(23) = {15,16};
Line(24) = {16,13};
Line(25) = {17,18};
Line(26) = {18,19};
Line(27) = {19,20};
Line(28) = {20,17};
Line(29) = {13,17};
Line(30) = {14,18};
Line(31) = {15,19};
Line(32) = {16,20};

//测线
Line(33) = {21,22};



//faces 围成面的线条有方向
Curve Loop(1) = {1, 2, 3, 4};
Curve Loop(2) = {5, 6, 7, 8};
Curve Loop(3) = {4, 13, -8, -16};
Curve Loop(4) = {14, 6, -15, -2};
Curve Loop(5) = {13, 5, -14, -1};
Curve Loop(6) = {16, -7, -15, 3};
Curve Loop(7) = {9, 10, 11, 12};
Curve Loop(8) = {17,-12, -20, 8};
Curve Loop(9) = {18, 10, -19, -6};
Curve Loop(10) = {17, 9, -18, -5};
Curve Loop(11) = {20, -11, -19, 7};

Curve Loop(12) = {21, 22, 23, 24};
Curve Loop(13) = {25, 26, 27, 28};
Curve Loop(14) = {29, -28, -32, 24};
Curve Loop(15) = {30, 26, -31, -22};
Curve Loop(16) = {29, 25, -30, -21};
Curve Loop(17) = {32, -27, -31, 23};

Plane Surface(1) = {1};
Plane Surface(2) = {2};
Plane Surface(3) = {3};
Plane Surface(4) = {4};
Plane Surface(5) = {5};
Plane Surface(6) = {6};
Plane Surface(7) = {7};
Plane Surface(8) = {8};
Plane Surface(9) = {9};
Plane Surface(10) = {10};
Plane Surface(11) = {11};
Plane Surface(12) = {12};
Plane Surface(13) = {13};
Plane Surface(14) = {14};
Plane Surface(15) = {15};
Plane Surface(16) = {16};
Plane Surface(17) = {17};

// —— 将外围六个面组成一个物理面组 "Outer" ——
// 顶面=1，底面=7，四个合并侧面=1001..1004
Physical Surface("Outer") = {1, 7, 3, 8, 4, 8, 5, 10, 6, 11};

//Vol
Surface Loop(1) = {1,6,3,5,4,2};
Surface Loop(2) = {2,11,8,10,9,7};
Surface Loop(3) = {12,17,14,16,15,13};

Volume(1) = {1};
Volume(2) = {2,3};
Volume(3) = {3};

//Add 3 physical tag
Physical Volume("Air",1) = {1};
Physical Volume("Ground",2) = {2};
Physical Volume("Abnormal",3) = {3};

MeshSize{PointsOf{Volume{:};}}=lc;
MeshSize{PointsOf{Volume{3};}}=lc/250;


Field[1] = Distance;
Field[1].CurvesList = {33};
Field[1].Sampling = 100;

Field[2] = Threshold;
Field[2].InField = 1;
Field[2].SizeMin = lc/300;
Field[2].SizeMax = lc;
Field[2].DistMin = 80;
Field[2].DistMax = 80000;

Background Field = 2;

Mesh.Algorithm = 5;

//Field[1] = Box;
//Field[1].Thickness = 2000;
//Field[1].VIn = lc/100;
//Field[1].VOut = lc;
//Field[1].XMin = 0;
//Field[1].XMax = 1000;
//Field[1].YMin = -500;
//Field[1].YMax = 500;
//Field[1].ZMin = 0;
//Field[1].ZMax = 1000;

//Field[2] = Box;
//Field[2].Thickness = 100;
//Field[2].VIn = lc/400;
//Field[2].VOut = lc;
//Field[2].XMin = -500;
//Field[2].XMax = 500;
//Field[2].YMin = -500;
//Field[2].YMax = 500;
//Field[2].ZMin = -50;
//Field[2].ZMax = 50;

//Field[3]  = Min;
//Field[3].FieldsList = {1,2};



//Mesh.MeshSizeExtendFromBoundary = 0;
//Mesh.MeshSizeFromPoints = 0;
//Mesh.MeshSizeFromCurvature = 0;

// This will prevent over-refinement due to small mesh sizes on the boundary.

// Finally, while the default "Frontal-Delaunay" 2D meshing algorithm
// (Mesh.Algorithm = 6) usually leads to the highest quality meshes, the
// "Delaunay" algorithm (Mesh.Algorithm = 5) will handle complex mesh size
// fields better - in particular size fields with large element size gradients:





